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ABSTRACT 


A  general-purpose  multivariable  functionalization  program  is  described 
in  terms  of  its  general  capabilities  and  the  information  which  it  affords 
the  user. 


foreword 


The  work  reported  herein  was  begun  under  WEPTASK  RT8801001/210-1/ 
54390001,  and  extended  under  WEPTASK  RM3773001/210-1/WS470B02 ,  in 
connection  with  the  development  of  a  general  purpose  functionalization 
program  for  se  in  ballistics. 
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INTRODUCTION 


The  problem  of  finding  an  analytical  formula  which  satisfactorily  repro¬ 
duces  experimental  data  or  data  which  are  the  end-product  of  more  or  less 
complex  computations  is  a  recurring  one;  and  the  number  of  necessary  or 
useful  functional  forms  is  boundless.  The  program  described  herein  affords 
the  user  a  wide  variety  of  functional  building  blocks  and  a  highly  flexi¬ 
ble  means  of  assembling  them  into  a  structure  which  will  serve  his  purpose. 
The  available  forms  are  listed  below. 


f 

M 

da 

hi 

db 

M 

dc 

1 .  a 

1 

2,  a 

coa  bxj 

cos  bxj 

-axj  sin  bxj 

3.  a 

aln  bxj 

sin  bxj 

axj  cos  bxj 

4 .  a 

e  bxi 

ebX1 

bxi 

axj  e 

/  5.  a 

b 

xi 

b 

X 

axj  log  xj 

6  a 

coan  bxj* 

cosn  bxi 

-na  xj  cos""1  bxj  *in  bxj 

7.  a 

sin"  bx j* 

8inn  bxj 

na  xj  aln"'1  bxj  cos  bxj 

8.  a 

In  bxj 

In  bxj 

a 

b 

9.  a 

coe  bxj  cos  cxj 

cos  bxj  cos  cxj 

-axj  sin  bxj  cos  cxj 

-axj  coa  bxj  sin  cx 

10.  a 

cos  bxj  aln  cxj 

cos  bxj  sin  cxj 

-axj  sin  bxj  sin  cxj 

axj  co*  bxj  co*  cx. 

11.  a 

bxi 

e  *  cos  cxj 

bx4 

e  1  cos  cxj 

bxj 

axj  e  1  co*  cxj 

bxi 

-axj  e  *in  cxj 

12  a 

b 

Xj  coa  cxj 

b 

Xj  CO*  CXj 

a(log  xj)  xj  co*  cxj 

-axj  (xj)  sin  cXj 

13.  a 

sin  bxj  sin  cx< 

sin  bxj  sin  cxj 

axj  cos  bxj  sin  cxj 

axj  *ln  bxj  cos  cxj 

14  a 

bxi 

e  sin  cxj 

bxj 

e  1  sin  cxj 

axj  ebx*  tin  cxj 

bxj 

axj  e  1  cos  cxj 

15  a 

xj  sin  cxj 

b  . 

Xj  sin  cxj 

a  (log  xj)  xb  sin  cxj 

axj  (xb)  co*  cxj 

18.  a 

eb*i  ec*j 

bx.  cxf 
e  1  e  J 

bxj  cx, 
a  Xj  e  e  J 

bx<  cxj 
axj  e  1  c  J 

17.  a 

eb*i  x^ 

bx  j  c 
e  xj 

c  b*i 
a  xj  Xj  e 

a  (log  Xj)  xj  ebxi 

18.  a 

b  c 

xi  Xj 

b  c 

Xi  x} 

a  xb  (log  Xj)  Xj 

a  xj  (log  xj)  Xj 

19.  a 

xj  Xj  cos  bx|,** 

20  a 

xj  xj  sin  bxj,** 

21  a 

xtbxJ 

a(xj  In  Xj)»jb*J 

*  n  s  2 ,  an  integer 

♦♦Available  for  use  In  the  linear  end  computet  lone!  mode*  only 


The  partial  derivatives  of  the  forms  with  respect  to  the  various  parameters 
of  fit  are  also  listed;  they  are  required  in  the  program  and  are  available, 
with  some  restrictions,  for  general  use.  The  letters,  a,  b,  £,  are  generic, 
of  course,  and  are  input  quantities  to  be  used,  determined  or  improved  as 
the  case  may  be.  The  letters  and  x^  are  input  variables,  dependent  or 
independent,  at  the  disposal  of  the  user.  The  number  of  input  variables, 
the  number  of  terms  in  the  fitting-function  and  the  number  of  data  points 
which  can  be  used  depend  on  the  dimensioning  of  the  program  to  fit  the 
internal  memory  of  the  computer. 

GENERAL  FLAN 


A  simple  way  to  grasp  the  general  capability  of  the  program  is  to  think 
of  it  in  terms  of 

1.  A  table  of,  say,  ten  columns  of  memory  of  which  any  or  all  can  be 
filled  with  input  variables  which  can,  so  far  as  the  machinery  of  the 
program  is  concerned,  be  used  as  dependent  or  independent  variables; 

2.  A  catalog  of  functions  and  their  derivatives  with  respect  to  the 
parameters  of  fit; 

3.  A  computational  "black  box"  which 

a.  Evaluates  a  selection  of  functions  from  the  catalog  and,  when 
it  is  appropriate,  their  derivatives,  at  the  values  of  selected  vari¬ 
ables  in  the  table  to  form  an  array  of  computed  values  which  will  be 
called  the  cp-matrix; 

b.  Multiplies  the  qj-matrix  by  its  transpose,  thus  forming  normal 
equations;  and 

c.  Solves  the  normal  equations  to  give  the  parameters  of  the  fit. 

Whenever  a  fit  is  made  the  printed  end-output  lists  the  parameters  which 
were  found,  the  average  of  the  absolute  values  of  the  residuals  of  the 
fit,  the  square  of  the  sum  of  the  residuals,  the  summed  squares  and 
variance  of  the  residuals;  the  average  of  the  absolute  values  of  the 
relative  error  at  each  point,  expressed  as  a  percent,  is  also  given 
with  the  numerically  largest  residual  and  the  point  at  which  it  occurs. 

The  phrase  "end-output"  is  used  above  because  there  are  several  optional 
print-outs;  those  just  listed  are  automatically  given  whenever  a  fit,  as 
distinguished  from  an  auxiliary  computation,  is  made. 

The  program  uses  the  central  machinery  for  producing  the  normal  matrix 
to  give  two  allied  modes  of  operation  for  functionalization:  a  linear . 
or  true  least-squares,  mode  and  a  non-linear  or  "differential  correction" 
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mode  which  will  herein  be  called  the  iterative  mode .  A  natural  by-product 
is  the  computational  mode  which  is  useful  both  as  an  auxiliary  in  pro¬ 
ducing  a  fit  and  producing  tables  of  functions  for  plotting  or  other  uses. 

THE  LINEAR  MODE 

The  linear  mode  is  :he  simplest  least-squares  type  of  operation  and  pro¬ 
duces  the  "least-squares  best"  fit  by  determining  the  coefficients  a4  in 
functions  of  the  form 

y  =  ^'Pi  +  +  ■■  ■  (i) 

in  which  the  cpt  can  be  drawn  in  any  manner  from  the  list  already  given. 

The  arguments  of  the  cp 1 s  may  be  any  input  variable.  It  will  be  noticed 
.that  the  listed  functions  individually  contain  either  one  or  two  variables, 
e.g.,  type  2  is  a  cos  bx^  in  which  the  subscript  i_  implies  that  £  can  be 
any  input  variable;  type  9  in  the  list  is  a  cos  bxt  cos  cx,  in  which  x^ 
and  Xj^  can  be  any  input  variables  or  the  same  variable  (see  type  6) .  The 
letters  a,  1),  and  £  are  arbitrary  input  quantities  in  any  case,  but  in  the 
linear  mode  the  "best"  a ' s  are  found  by  the  program.  The  customer  decides 
in  which  column  of  the  data-table  he  wishes  to  place  a  particular  variable 
and  then  calls  for  it  by  number. 

A  series  of  cosines  of  a  single  argument,  for  example,  would  be  set  up  by 
repetitive  use  of  type  2  in  the  function-list.  Suppose  that  the  values  of 
the  common  argument  of  the  trigonometric  terms  have  been  put  into  column 
five,  then  the  cosine  series 

y  =  cos  x  +  cos  2x  +  cos  3x  +  cos  7x  (2) 

would  be  called  from  the  program  by  a  card  showing 

2  5,  2  5,  2  5,  2  5 

based  on  type  2,  i.e.,  a  cos  bx,  ,  and  variable  number  5.  The  desired 
harmonics  would  be  called  by  another  card  showing 

1  1,  1  2,  1  3,  1  7 

which  give  the  requisite  values  of  a  and  b  in  a  cos  bt .  Although  the 

inputs  a_y  b^,  c_  in  the  forms  are  arbitrary,  some  values  must  always  be 

given. 

Examination  oi  the  list  of  functions  shows  that  single  and  double  Fourier 
series  and  the  redoubtable  least -squares  polynomial  are  immediately  avail¬ 
able.  Certain  options  to  be  discussed  later  allow  approximations  in  terms 
of  various  types  of  polynomials  such  as  those  of  Hermite,  Legendre, 
Chebychev,  etc.  The  number  of  terms  of  such  expansions  depends  primarily 
on  the  number  oi  slots  allowed  for  input  variables.  Various  standard 
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multivariable  expansions  can  be  concocted.  The  coding  of  the  program  has 
been  laid  out  in  such  a  way  as  to  allow  the  list  of  basic  functions  to  be 
altered  on  an  ad  hoc  basis.  For  example,  the  general-purpose  program  was 
readily  adapted  to  allow  175  terms  of  double  Fourier  type,  any  or  all  of 
which  might  be  modified  by  a  positive  or  negative  power  of  a  non-angular 
variable;  this  adaptation  handles  a  maximum  of  9000  data-points. 

Examples  of  types  of  operation  which  can  be  followed  will  be  given  later. 

THE  ITERATIVE  MODE 


The  iterative  mode  of  the  program  is  a  succession  of  applications  of  the 
differential  correction  technique  allowing  the  determination  not  only  of 
the  linear  coefficient,  a,  appearing  in  the  list  of  forms,  but  the  para¬ 
meters  1j  and  £  as  well. 

The  technique  used  (Reference  1)  calls  for  the  first  differential  of  the 
fitting-function  with  respect  to  the  various  parameters;  and  these  are 
catalogued  internally,  so  to  speak,  in  the  program:  when  the  control 
calling  for  the  iterative  mode  is  used  the  program  automatically  supplies 
the  derivatives  which  are  required. 

The  quantities  which  are  found  by  solving  the  normal  equations  are  used 
as  corrections  and  are  applied  to  the  input  values  of  the  parameters; 
a  new  set  of  normal  equations  can  be  formed  using  the  altered  values  to 
find  new  increments,  and  so  on.  The  process  can  be  terminated  after  a 
preassigned  tolerance  is  met  or  after  a  preassigned  number  of  iterations 
if  the  tolerance  is  not  met.  The  initial  values  of  the  parameters  have 
co  be  reasonably  good,  however,  or  the  function  is  likely  to  blow  up. 

A  few  ways  of  getting  the  process  started  on  a  reasonably  sound  basis 
will  be  discussed  in  a  later  section. 

THE  COMPUTATIONAL  MODE 


The  computational  mode  involves  no  normal  equations  or  solutions  thereof. 
It  is  intended  to  extend  the  number  of  functions  which  are  available  to 
the  user  on  the  one  hand,  or  to  generate  funct ion- values  to  be  used  for 
purposes  other  than  fitting.  Generally  the  computational  mode  allows  one 
to  list  or  punch  on  cards  practically  any  set  of  numbers  which  could 
enter  into  the  computation  of  a  fit. 

OPTIONS  AND  OUTPUT 

It  is  helpful  in  discussing  options  and  output  to  speak  of  input  "columns" 
ai  i  the  columns  of  the  p-matrix,  i.e.,  the  columns  of  the  matrix  of  the 
equations  of  condition. 
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I .  Operational  Options 

A.  Linear  least  scares 

B.  Differential  correction,  iterative  mode 

C.  "Sum  cp";  this  option  tells  the  computer  to  add  the  elements  of 
the  rows  of  the  <p-matrix  row  by  row  for  future  use,  without  completing 
the  leas t -squares  routine. 


II •  Print-Only  Options 
A. 


Long 

print ; 

1.  a 

.  The 

b 

The 

r 

.  The 

d 

.  The 

e 

.  The 

f 

.  The 

and 

g 

The 

absolute  values  of  (residual/input  value)  x  100. 

(Note:  This  information  is  given  whenever  a  fit  is  made.) 

a.  The  values  of  the  function  be^ng  fitted, 

b.  The  values  of  the  input  variables  used  in  the  fit, 

c .  The  residuals . 

The  cp-matrix. 

The  matrices  of  the  normal  equations 

(Note:  This  information  is  given  after  each  step  in  the 
iterative  mode,  with  the  "correction"  to  each  parameter, 
as  well  as  the  corrected  parameters.) 
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The 

normal  print;  this 

gives 

1. 

Item  A . 1 . 

above 

2. 

Item  A. 2. 

above 

The 

short  print;  this 

gives 

1  . 

Item  A . 1 . 

above . 

(The  input  controls  are  always  printed.) 

HI  •  Print-Punch-Store  Option 

The  following  items  may  be  listed,  card-punched  or  stored  as  input 
variables  for  subsequent  use: 

A.  The  fitted  function;  i.e.,  the  computed  values  of  the  fit; 

B.  The  row -by-row  sum  of  the  elements  of  the  cp-matrix; 

C.  Any  column  of  the  cp-matrix. 

IV.  Miscellaneous  Options 

A.  Segmentation  Option.  This  option  allows  the  user  to  select  any 
block  of  consecutive  data  points  for  fitting,  e.g.,  points  1-25,  or  31-268. 

B.  Point-Selection  Option.  This  option  allows  every  n1-*1  point  of  the 
entire  set  of  input  data  or  segment  thereof  to  be  used  in  a  fit;  the  resid¬ 
uals,  if  called  for,  are  computed  over  the  entire  block. 

C.  Carry-Over  Options.  These  options  allow: 

1.  A  particular  functional  form  to  be  used  over  several  segments 
of  data  without  a  full  set  of  controls  for  each  segment; 

2.  The  parameters  as  determined  from  one  fit  to  be  used  (in  the 
same  functional  form)  in  the  next  operation. 

EXAMPLES  AND  DISCUSSION 


The  examples  of  uses  to  which  the  program  has  been  put  will  not  be  given 
in  terms  of  input  or  output  numbers,  but  in  terms  of  the  ways  in  which 
the  problems  were  handled. 

In  one  case  a  shotgun  approach  to  the  following  problem  was  used.  Given 
some  90  sets  of  72  similar  but  different  physical  measurements  all  made 
at  equal  intervals  of  a  single  independent  angular  variable  it  was  desired 
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to  investigate  whether  these  sets  could  be  fitted  with  a  Fourier  type  of 
function  of  sines  and  cosines,  allowing  as  many  as  four  harmonics.  The 
version  of  the  program  which  was  used  allowed  a  maximum  of  700  data  points 
in  each  of  ten  allowable  variable-columns  (and,  incidentally,  a  maximum 
of  30  terms  in  the  linear  mode). 

Seventy-two  of  the  sets  were  handled  in  one  pass  through  the  computer  by 
loading  each  of  9  input  columns  with  9  complete  sets  of  measured  values, 
and  column  10  with  9  identical,  complete  sets  of  values  of  the  angular 
variable;  using  the  segmentation  option  and  a  carry-over  option,  all  of 
the  81  stored  sets  were  fitted  successively  with  first  harmonic,  first 
and  second  harmonics,  etc.  The  resulting  fits  were  card-punched  for  later 
use  in  a  plotting  routine. 

A  simple  example  using  the  iterative  mode  is  furnished  by  the  problem  of 
functionalizing  the  speed  of  a  projectile  in  terms  of  its  distance  from 
the  muzzle  of  the  gun.  The  function  used  was 


v  being  speed  and  s  distance  from  muzzle.  It  was  desired  in  this  case 
to  find  suitable  values  of  a  and  k.  From  the  general  knowledge  avail¬ 
able,  a  good  value  of  k  was  at  hand  and  a  plausible  value  of  v  could  be 
estimated  from  the  data,  so  that  the  sometimes  knotty  problem  of  finding 
initial  values  for  the  iterative  process  was  easily  solved. 


A  problem  which  allows  the  program  a  little  more  scope  is  furnished  by 
a  particular  approach  to  the  functionalization  of  bombing-table  ranges. 

The  ranges  to  be  reproduced  are  taken  as  functions  of  height,  speed  and 
angle  of  dive  or  climb  of  the  aircraft  at  the  time  of  release  of  the  bomb. 
A  function  which  serves  nicely  for  one  type  of  bomb  involves  finding 
values  of  a1  ,  Og  ,  .  .  .  ,  ag  in 


R  =  a1  v  cos  ag  0  exp(a3  +  a4v  +  ag 


A  step-by-step  technique  for  finding  the  six  parameters  in  this  function 
will  be  given;  the  reasons  for  the  operations  are  not  relevant.  Using 
data  at  9  =  0  (level  flight), 


1  Compute  R/V  (type  18  in  the  list  of  functions)  and  store; 

2.  Compute  In  R/V  and  store  (type  8); 

3.  Fit  In  R/V  -  ,t,+a4+as^-i-+a6^:-  (types  1,  5,  18); 

Using  data  based  on,  say,  dive  angles,  and  the  parameters  found  in 
Step  3; 
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4.  Compute  exp(a3  +a4v  +  .it;  j —  +  ae  ~  )  and  store; 


5.  Compute  R/’V  exp(a3  +  a4v  +  —  +  ag  —  )]  and  store; 

h  v 


6.  Fit  the  result  of  step  5  with  aj  cos  as0  (type  2)  using  the  iterative 
mode  to  evaluate  and  a^  . 

If  the  total  number  of  data  points  does  not  exceed  the  capacity  of  the 
program,  the  process  can  be  completed  in  a  single  pass  through  the  computer. 

It  will  be  apparent  that  so  far  as  the  linear  mode  of  operation  is  con¬ 
cerned,  none  of  the  products,  types  6,  7  and  9-20,  is  logically  necessary; 
but  in  order  to  keep  the  entire  program  with  its  many  arrays  within  the 
internal  memory  of  the  computer,  and  at  the  same  time  amendable  to  FORTRAN 
coding,  one  has  to  start  and  stop  somewhere.  If  ten  columns  of  n  locations 
are  allowed  for  input  or  computed  quantities  and  a  maximum  of  10  terms  in 
the  fit,  41n  locations  are  gone  without  any  coding.  However,  there  will 
be  available  for  storing  composite  functions  ten  columns  less  the  number 
of  input  variables. 


The  picture  is  somewhat  different  in  the  case  of  the  iterative  mode. 

Each  parameter  of  fit  gives  rise  to  a  column  of  the  rp-array  by  way  of 
the  differential.  Hence  30  available  columns  in  the  cp-array  permit  the 
finding  of  thirty  parameters  which  may  be  involved  in  ten  three-parameter 
functionals,  or  fifteen  two-parameter  types,  or  however.  The  whole  scheme 
could  doubtless  be  greatly  improved  at  the  level  of  computer-system  de¬ 
velopment  . 


Use  of  the  iterative  mode  must  be  used  with  caution.  The  technique  of 
differential  correction  calls  for  the  evaluation  of  the  derivatives, 
shown  in  the  list  of  functions,  at  some  set  of  values  of  the  parameters 
at  each  step  of  the  iteration  as  they  are  inexorably  ground  out  by  the 
computer;  and  the  only  control  one  has  over  the  process  is  in  the  choice 
of  the  first  set,  which  is  an  input,  and  this  choice  is  accordingly  very 
important.  One-term  fits  which  involve  powers  or  exponentials  of  a 
variable  or  products  thereof  can  be  managed  by  selecting  a  data  point 
for  each  parameter  and  setting  up  two  or  three  linear  equations  by 
transforming  to  logarithms;  these  equations  can  then  be  solved  for  the 
parameters.  Lacking  any  clue  at  all,  a  cut  and  try  approach  can  be 
taken  by  selecting  a  data  point  for  each  parameter  -  or  several  such 
sets  if  there  is  a  large  amount  of  data  covering  wide  ranges  of  values 
of  the  variables  -  and  producing  a  set  of  values  for  trial  or  examinati'  i 
by  callinp  for  a  single  iteration  based  on  input  parameters  set  equal 
for  unity. 

A  systemmatic  analysis  of  the  questions  involved  here  does  not  seem  tc 
have  been  carried  out  ;  it  would  be  a  major  enterprise,  but  would  be  a 
valuable  advancement  of  the  art. 
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